#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <string.h>
#include <mpi.h>
#include <limits.h>
#include <sys/cdefs.h>
#include <tuple>
#include <type_traits>
#include <unistd.h>
#include <libgen.h>
#include <signal.h>
#include <setjmp.h>
#include <execinfo.h>
#include "charmm_param.h"
#include "esmd_types.h"
#include "htab_cpp.hpp"
#include "linux/limits.h"
#include "log.h"
#include "io_pdb.h"
#include "io_psf.h"
#include "io_charmm_inp.h"
#include "elemtab.h"
#include "cell.h"
#include "dimension.h"
#include "bonded.h"
#include "comm.h"
#include "nonbonded.h"
#include "coul_msm.h"
#include "esmd_conf.h"
#include "charmm_force.h"
#include "force_field.h"
#include "param_mapper.hpp"
#include "unitconv.hpp"
#include "velocity.h"
#include "verlet.h"
#include "shake.h"
#include "rigid.h"
//#include "rattle.h"
#include "timer.h"
#include "minimize.h"
#include "nose_hoover.h"
#include "bt.h"
#include "memory_cpp.hpp"
#undef min
#undef max
#include "path.hpp"
#include "rigid.hpp"
#include "lincs.hpp"
#include "shake.hpp"
#include "trajectory.hpp"
#include "listed.hpp"
#include "gmxtop.hpp"
#include "spring_self.h"
#include "spring_com.h"
#ifdef __sw__
#include "sw/swarch.h"
#include "qthread.h"
#endif
extern void compute_listed_forces(cellgrid_t *grid, void *param, mdstat_t *stat);
size_t file_size(FILE *f) {
  fseek(f, 0, SEEK_END);
  size_t off = ftell(f);
  fseek(f, 0, SEEK_SET);
  return off;
}
extern void cell_bfs(cellgrid_t *grid);
extern void part_cellgrid(cellgrid_t *grid, int nn, real rcut, real skin, const box<real> &gbox, const box<real> &lbox);
extern void forward_comm_bond(cellgrid_t *grid, mpp_t *mpp);
static jmp_buf time_machine;
void supress_segv(int sig) {
  puts("backtrace failed!");
  longjmp(time_machine, 1);
}
int backtrace_safe(void **buf, int n) {
  struct sigaction sa;
  sigaction(SIGSEGV, NULL, &sa);
  __sighandler_t old_handler = sa.sa_handler;
  sa.sa_handler = supress_segv;
  sigaction(SIGSEGV, &sa, NULL);
  int error = setjmp(time_machine);
  int ntrace;
  if (!error) {
    ntrace = backtrace(buf, n);
    puts("backtrace success");
  } else {
    ntrace = 0;
  }
  buf[ntrace] = 0;
  sa.sa_handler = old_handler;
  sigaction(SIGSEGV, &sa, NULL);
  return ntrace;
}
int gr;

#ifdef __sw__
#include <athread.h>
#endif
void sigsegv(int sig) {
  printf("sig=%d\n", sig);
  void *buf[100];
  int ntrace = backtrace(buf, 100);
  puts("after backtrace");
  fflush(stdout);

  printf("%d-%d %d\n", gr, sig, ntrace);
  fflush(stdout);

  for (int i = 0; i < ntrace; i++) {
    printf("%p\n", buf[i]);
  }
#ifdef __sw_has_speio__
  for (int i = 0; i < 64; i++) {
    long pc = SPEIO_READ_REG(0, i, SPE_IBXIO_GPC);
    printf("%02d: %p\n", i, pc);
  }
#endif
  fflush(stdout);
  exit(1);
  // while (1)
  //   ;
}
void sigbt(int sig) {
  char fn[4096];
  sprintf(fn, "backtrace.%06d.txt", gr);
  FILE *f = fopen(fn, "w");
  void *buf[100];
  int ntrace = backtrace_safe(buf, 100);
  for (int i = 0; i < ntrace; i++) {
    fprintf(f, "%p\n", buf[i]);
  }
  fclose(f);
}
__always_inline long displace_tag(long tag, long off, long total) {
  long ret = tag + off;
  while (ret < 0)
    ret += total;
  while (ret >= total)
    ret -= total;
  return ret;
}

long init_base_cell_namd(cellgrid_t *grid, elemtab_t *elemtab, mpp_t *mpp, esmd_config_t *conf) {
  char pdbpath[2048], psfpath[2048];
  cat_expand_path(pdbpath, conf->root, conf->coordinates);
  cat_expand_path(psfpath, conf->root, conf->structure);
  /*get x from pdb*/
  
  long natoms;
  double(*x)[3];
  if (mpp->pid == 0) {
    pdb_data_t *pdb = read_pdb(fopen(pdbpath, "r"));
    natoms = pdb->natom;
    x = esmd::allocate<double[3]>(natoms, "main/x");
    for (int i = 0; i < natoms; i++) {
      x[i][0] = pdb->atom[i].x;
      x[i][1] = pdb->atom[i].y;
      x[i][2] = pdb->atom[i].z;
    }
    free_pdb(pdb);
  }
  MPI_Bcast(&natoms, 1, MPI_LONG_LONG, 0, MPI_COMM_WORLD);
  if (mpp->pid != 0) {
    x = esmd::allocate<double[3]>(natoms, "main/x");
  }
  MPI_Bcast(x, natoms * 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  //free_pdb(pdb);

  int *t = esmd::allocate<int>(natoms,"main/type");
  real *q = esmd::allocate<real>(natoms,"main/charge");
  real *mass = esmd::allocate<real>(natoms,"main/mass");

  /*q and t from psf*/
  psf_data_t *psf;
  if (mpp->pid == 0) {
    psf = read_psf(fopen(psfpath, "r"));
  } else {
    psf = esmd::allocate<psf_data_t>(1, "psf");
  }
  MPI_Bcast(&psf->natom, 8, MPI_LONG_LONG, 0, MPI_COMM_WORLD);
  if (mpp->pid > 0) {
    //printf("%d %d\n", psf->nbond, psf->nimpr);
    psf->atom = esmd::allocate<psf_atom_t>(natoms, "psf/atoms");
    psf->bond = esmd::allocate<bond_t>(psf->nbond, "psf/bonds");
    psf->angle = esmd::allocate<angle_t>(psf->nangle, "psf/angles");
    psf->tori = esmd::allocate<dihed_t>(psf->ntori, "psf/toris");
    psf->impr = esmd::allocate<dihed_t>(psf->nimpr, "psf/imprs");
  }
  MPI_Bcast(psf->atom, sizeof(psf_atom_t) * natoms, MPI_BYTE, 0, MPI_COMM_WORLD);
  MPI_Bcast(psf->bond, psf->nbond * 2, MPI_LONG_LONG, 0, MPI_COMM_WORLD);
  MPI_Bcast(psf->angle, psf->nangle * 3, MPI_LONG_LONG, 0, MPI_COMM_WORLD);
  MPI_Bcast(psf->tori, psf->ntori * 4, MPI_LONG_LONG, 0, MPI_COMM_WORLD);
  MPI_Bcast(psf->impr, psf->nimpr * 4, MPI_LONG_LONG, 0, MPI_COMM_WORLD);
  //puts("after bcast");
  for (int i = 0; i < natoms; i++) {
    int ocnt = elemtab->cnt;
    t[i] = elemtab_map(elemtab, psf->atom[i].type);
    if (mpp->pid == 0 && elemtab->cnt != ocnt) {
      printf("%d %s %f\n", t[i], psf->atom[i].type, psf->atom[i].mass);
    }
    q[i] = psf->atom[i].charge;
    mass[i] = psf->atom[i].mass;
  }
  vec<real> glen = conf->cell.hi - conf->cell.lo;
  // vecsubv(glen, conf.cell.hi, conf.cell.lo);
  if (mpp->pid == 0) {
    printf("Building bond graph\n");
  }
  // if (mpp->pid == 0) printf("maxdb=%f %d %d\n", maxdb, mi, mj);
  bond_graph_t graph;
  build_graph(&graph, psf);
  esmd::deallocate(psf->bond);
  esmd::deallocate(psf->angle);
  esmd::deallocate(psf->tori);
  impr_index_t imidx;
  index_impr(&imidx, &graph, psf);
  //free_psf(psf);
  if (mpp->pid == 0) {
    printf("Bcasting params\n");
  }

  // long(*excls)[MAX_EXCLS_ATOM] = esmd::allocate<long[MAX_EXCLS_ATOM]>(natoms, "main/excls");
  // long(*scales)[MAX_SCALS_ATOM] = esmd::allocate<long[MAX_SCALS_ATOM]>(natoms, "main/scals");
  // long(*chain2)[MAX_SCALS_ATOM][2] = esmd::allocate<long[MAX_SCALS_ATOM][2]>(natoms, "main/chain2");
  //bfs_excls_chain2(natoms, excls, chain2, scales, &graph);

  if (mpp->pid == 0) {

    printf("Building cells\n");
  }

  build_cells_bondonly(grid, 2, conf->rcut, conf->skin, &mpp->gbox, &mpp->lbox, natoms, x, q, t, mass, &graph, &imidx, OVF_SKIP);
  // }
  grid->gbox = mpp->gbox;
  // boxcpy(grid.gbox, gbox);
  if (mpp->pid == 0) {
    printf("Creating velocity\n");
  }
  
  create_velocity(grid, conf->temperature, mpp);
  return 0;
}
template<long (FALLBACK_INIT)(cellgrid_t *grid, elemtab_t *elemtab, mpp_t *mpp, esmd_config_t *conf)>
long init_from_restart(cellgrid_t *grid, elemtab_t *elemtab, mpp_t *mpp, esmd_config_t *conf) {
  if (conf->restart != NULL) {
    part_cellgrid(grid, 2, conf->rcut, conf->skin, mpp->gbox, mpp->lbox);
    char restart_path[PATH_MAX];
    expand_path(restart_path, conf->root, "%s", conf->restart);
    long istep = read_restart(restart_path, mpp, grid, elemtab);
    printf("restarting from step %ld\n", istep);
    if (istep >= 0)
      return istep;
  }
  return FALLBACK_INIT(grid, elemtab, mpp, conf);
}
void init_charmm_param(charmm_param_t *param, mpp_t *mpp, elemtab_t *elemtab, esmd_config_t *conf){
  /*process inp*/
  charmm_rawparam_t *inp;
  char inppath[PATH_MAX];
  cat_expand_path(inppath, conf->root, conf->parameters);
  if (mpp->pid == 0) {
    inp = read_inp(fopen(inppath, "r"));
  } else {
    inp = esmd::allocate<charmm_rawparam_t>(1, "main/inp");
  }
  MPI_Bcast(&inp->nbond, 5, MPI_INT, 0, MPI_COMM_WORLD);
  if (mpp->pid > 0) {
    inp->bond = esmd::allocate<charmm_raw_bond_param_t>(inp->nbond, "inp/bond");
    inp->angle = esmd::allocate<charmm_raw_angle_param_t>(inp->nangle, "inp/angle");
    inp->tori = esmd::allocate<charmm_raw_dihed_param_t>(inp->ntori, "inp/tori");
    inp->impr = esmd::allocate<charmm_raw_dihed_param_t>(inp->nimpr, "inp/impr");
    inp->nonb = esmd::allocate<charmm_raw_nonb_param_t>(inp->nnonb, "inp/nonb");
  }
  MPI_Bcast(inp->bond, sizeof(charmm_raw_bond_param_t) * inp->nbond, MPI_BYTE, 0, MPI_COMM_WORLD);
  MPI_Bcast(inp->angle, sizeof(charmm_raw_angle_param_t) * inp->nangle, MPI_BYTE, 0, MPI_COMM_WORLD);
  MPI_Bcast(inp->tori, sizeof(charmm_raw_dihed_param_t) * inp->ntori, MPI_BYTE, 0, MPI_COMM_WORLD);
  MPI_Bcast(inp->impr, sizeof(charmm_raw_dihed_param_t) * inp->nimpr, MPI_BYTE, 0, MPI_COMM_WORLD);
  MPI_Bcast(inp->nonb, sizeof(charmm_raw_nonb_param_t) * inp->nnonb, MPI_BYTE, 0, MPI_COMM_WORLD);
  if (mpp->pid == 0) {
    printf("Making params\n");
  }

  build_param(param, inp, elemtab, conf->rcut, conf->rswitch, conf->coulconst, conf->coulscaling);
  param->nonb.coultype = conf->coultype;
}

void init_signals(){
  signal(SIGABRT, sigsegv);
  struct sigaction sa;
  sigaction(SIGSEGV, NULL, &sa);
  __sighandler_t old_handler = sa.sa_handler;
  sa.sa_handler = sigsegv;
  sa.sa_flags |= SA_NODEFER;
  sigaction(SIGSEGV, &sa, NULL);
  sigaction(SIGILL, &sa, NULL);
  signal(SIGUSR1, sigbt);
}

__always_inline std::tuple<tagint, vec<real>*> read_coordinates(mpp_t &mpp, const esmd_config_t &conf){
  char pdbpath[PATH_MAX];
  cat_expand_path(pdbpath, conf.root, conf.coordinates);
  vec<real> *x;
  tagint natom;
  if (mpp.pid == 0){
    pdb_data_t *pdb = read_pdb(fopen(pdbpath, "r"));
    natom = pdb->natom;
    x = esmd::allocate<vec<real>>(natom, "main/x");
    for (tagint i = 0; i < natom; i ++){
      x[i].x = pdb->atom[i].x;
      x[i].y = pdb->atom[i].y;
      x[i].z = pdb->atom[i].z;
    }
    free_pdb(pdb);
  }
  natom = mpp.bcast_var(natom);

  if (mpp.pid != 0){
    x = esmd::allocate<vec<real>>(natom, "main/x");
  }
  mpp.bcast(x, natom);
  return std::make_tuple(natom, x);
}

void read_top(gmxtop &top, mpp_t &mpp, const esmd_config_t &conf){
  char toppath[2048];
  cat_expand_path(toppath, conf.root, conf.structure);
  FILE *topfile = fopen(toppath, "r");
  if (mpp.pid == 0){
    top.verbose = true;
  }
  top.parse_file(topfile);
  top.convert_settle();
}

void init_atoms_in_grid(cellgrid_t &grid, int nxatom, vec<real> *x, gmxtop &top){
  // grid.tag_map = esmd::create<esmd::htab<tag_key, vec<int>>>(1, "tag_map", "tag_map");
  // grid.type_map = esmd::create<esmd::htab<tag_key, int>>(1, "type_map", "type_map");
  long mol_cur = 0;
  grid.tag_map = new esmd::htab<tag_key, vec<int>>("tag_map");
  grid.type_map = new esmd::htab<tag_key, int>("type_map");
  auto &type_map = *grid.type_map;
  auto &tag_map = *grid.tag_map;
  tagint offset = 0;
  auto glen = grid.gbox.hi - grid.gbox.lo;
  for (auto &mol : top.molecules){
    auto &moltype = top.get_moltype(mol.name);
    for (int imol = 0; imol < mol.cnt; imol ++){
      for (int i = 0; i < moltype.atoms.size(); i ++){
        gmxtop::atom &atom = moltype.atoms[i];
        int itype = type_map.setdefault(atom.type.ival, type_map.size());

        vec<real> xi = x[offset + i];
        while (xi.x <  grid.gbox.lo.x) xi.x += glen.x;
        while (xi.y <  grid.gbox.lo.y) xi.y += glen.y;
        while (xi.z <  grid.gbox.lo.z) xi.z += glen.z;
        while (xi.x >= grid.gbox.hi.x) xi.x -= glen.x;
        while (xi.y >= grid.gbox.hi.y) xi.y -= glen.y;
        while (xi.z >= grid.gbox.hi.z) xi.z -= glen.z;
        if (!grid.lbox.contains(xi)) continue;
        
        vec<real> cf = ((xi - grid.lbox.lo) * grid.rlen);
        vec<int> c = {(int)floor(cf.x), (int)floor(cf.y), (int)floor(cf.z)};

        celldata_t &cell = grid.cell_at(c.x, c.y, c.z);


        cell.x[cell.natom] = xi - cell.basis;
        cell.xinit[cell.natom] = xi;
        cell.q[cell.natom] = atom.q;
        cell.mass[cell.natom] = atom.mass;

        cell.t[cell.natom] = itype;
        cell.rmass[cell.natom] = 1 / atom.mass;
        cell.tag[cell.natom] = offset + i;
        cell.mol[cell.natom] = mol_cur;
        cell.natom ++;
        tag_map[offset + i] = c;
      }
      offset += moltype.atoms.size();
      mol_cur += 1;
    }
  }
}
#include <typeinfo>
template<typename T, int IFUNC_SELS>
void gmx_to_topogrid(cellgrid_t &grid, gmxtop &top, mpp_t &mpp){
  listed_force<T> list;
  param_mapper<T> param_mapper("param mapper");
  top.export_listed(list, param_mapper, IFUNC_SELS);
  //if (mpp.pid == 0) printf("#listed_force<%s>=%ld\n", typeid(T).name(), list.topo.size());
  grid.topo->allocate_and_distribute(grid, *grid.tag_map, list);
  if (mpp.pid == 0) printf("#listed_force<%s>=%ld\n", typeid(T).name(), list.topo.size());
  
}
void init_topology(cellgrid_t &grid, gmxtop &top, mpp_t &mpp, esmd_config_t &conf){
  top.build_exclusions();

  top.unitconv(units::Ang, units::kcal_mole, units::rad);
  top.sort_rigids();
  
  listed_force<rigid_param> rigids;
  param_mapper<rigid_param> rigids_mapper("rigids mapper");
  if (conf.rigidtype != RIGID_NONE){
    top.find_clusters(rigids, rigids_mapper);
  }
  if (mpp.pid == 0) printf("number of rigids=%d %d\n", rigids.topo.size(), rigids.nparam);
  gmx_to_topogrid<harmonic_bond_param, 1<<1>(grid, top, mpp);
  gmx_to_topogrid<harmonic_angle_param, 1<<1>(grid, top, mpp);
  gmx_to_topogrid<urey_bradly_param, 1<<5>(grid, top, mpp);
  gmx_to_topogrid<cosine_torsion_param, (1<<1|1<<9)>(grid, top, mpp);
  gmx_to_topogrid<ryckaert_bellmans_param, 1<<3>(grid, top, mpp);
  gmx_to_topogrid<harmonic_improper_param, 1<<2>(grid, top, mpp);
  // listed_force<harmonic_bond_param> bonds;
  // param_mapper<harmonic_bond_param> bonds_mapper("bonds mapper");
  // top.export_listed(bonds, bonds_mapper, 1<<1);
  // if (mpp.pid == 0) printf("#bonds=%lu\n", bonds.topo.size());
  // listed_force<urey_bradly_param> angles;
  // param_mapper<urey_bradly_param> angles_mapper("angles mapper");
  // top.export_listed(angles, angles_mapper, 1<<5);
  // if (mpp.pid == 0) printf("#angles=%lu\n", angles.topo.size());
  // listed_force<cosine_torsion_param> tors;
  // param_mapper<cosine_torsion_param> tors_mapper("tors mapper");
  // top.export_listed(tors, tors_mapper, 1<<9);
  // if (mpp.pid == 0) printf("#tors=%lu\n", tors.topo.size());

  // listed_force<harmonic_improper_param> imprs;
  // param_mapper<harmonic_improper_param> imprs_mapper("imprs mapper");
  // top.export_listed(imprs, imprs_mapper, 1<<2);
  // if (mpp.pid == 0) printf("#imprs=%lu\n", imprs.topo.size());

  listed_force<special_pair> special_pairs;
  top.export_special_pairs(special_pairs);

  // if (mpp.pid == 0) printf("Total bonded param size: %ld\n", bonds.nparam * sizeof(*bonds.param) + angles.nparam * sizeof(*angles.param) + tors.nparam * sizeof(*tors.param) + imprs.nparam * sizeof(*imprs.param));
  //for (auto &atomtype : top.atom_types){
  //  printf("%s %f %f\n", atomtype.name.str, atomtype.epsilon, atomtype.sigma);
  //}
  grid.topo->allocate_and_distribute(grid, *grid.tag_map, special_pairs);
  if (conf.rigidtype != RIGID_NONE) {
    grid.topo->allocate_and_distribute(grid, *grid.tag_map, rigids);
  }
}

void init_params(charmm_param_t &param, esmd::htab<tag_key, int> &type_map, gmxtop &top, esmd_config_t &conf){
  esmd::htab<tag_key, std::tuple<real, real, real, real>> nbp_map("nonbonded params");
  top.create_vdw14_params(nbp_map);
  esmd::list<int> types("gmxtop types");
  param.nonb.vdw_param = esmd::allocate<vdw_param_t>(type_map.size(), "nonbonded params");
  param.ntypes = type_map.size();
  auto &nonb = param.nonb;
  std::tie(nonb.ntypes, nonb.coultype, nonb.coul_const, nonb.order, nonb.rcut, nonb.rsw, nonb.scale14) = std::tie(param.ntypes, conf.coultype, conf.coulconst, conf.msmorder, conf.rcut, conf.rswitch, conf.coulscaling);
  auto factor_sig = (cbrt(sqrt(2.0)) * units::nm / units::Ang).value;
  auto factor_eps = (units::kJ / units::mole / units::kcal_mole).value;
  
  for (auto &ent : type_map) {
    int it = ent.val;
    long gt = ent.key;
    auto &pm = param.nonb.vdw_param[it];
    elemtag t;
    t.ival = gt;

    //printf("%f %f %f %f\n", std::get<0>(nbp_map[gt]), std::get<1>(nbp_map[gt]), std::get<2>(nbp_map[gt]), std::get<3>(nbp_map[gt]));
    std::tie(pm.eps, pm.sig, pm.eps14, pm.sig14) = nbp_map[gt];
    pm.sig *= factor_sig * 0.5;
    pm.sig14 *= factor_sig * 0.5;
    pm.eps *= factor_eps;
    pm.eps14 *= factor_eps;
    pm.eps = sqrt(pm.eps);
    pm.eps14 = sqrt(pm.eps14);

  }
  //esmd::htab<tag_key, int> type_map("type map");

}

int main(int argc, char **argv) {
  // char *tmp[9];
  // gets(tmp);
  MPI_Init(&argc, &argv);
  int rank, nproc;
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Comm_size(MPI_COMM_WORLD, &nproc);
  gr = rank;
  init_prof_(&rank);
  start_prof();
  timer_init();
#ifdef __sw__
  qthread_init();
#endif
  init_signals();

  /*read config*/
  if (rank == 0) puts("reading config");
  esmd_config_t conf;
  parse_esmd_config_mpi(&conf, argv[1], rank);

  /*initialize global path envvars*/
  setenv("WORKDIR", getenv("PWD"), 1);
  setenv("CONFDIR", dirname(argv[1]), 1);
  assert(conf.structure != NULL && conf.coordinates != NULL);

  /*init box*/
  if (rank == 0) puts("initializing:");
  if (rank == 0) puts("  box");
  box<real> gbox;
  gbox = conf.cell;

  /*init communicator*/
  if (rank == 0) puts("  comm");
  mpp_t mpp;
  comm_init(&mpp, MPI_COMM_WORLD, &gbox);

  /*read coordinates*/
  if (rank == 0) puts("  coords");
  vec<real> *x;
  tagint natoms;
  std::tie(natoms, x) = read_coordinates(mpp, conf);
  /*read top*/
  if (rank == 0) puts("  topology");
  gmxtop top;
  read_top(top, mpp, conf);
  /*grid*/
  if (rank == 0) puts("  grid");
  cellgrid_t grid;
  part_cellgrid(&grid, 2, conf.rcut, conf.skin, mpp.gbox, mpp.lbox);
  init_atoms_in_grid(grid, natoms, x, top);
  /*topology*/
  if (rank == 0) puts("  topology grid");
  topology_grids topo_grid;
  grid.topo = &topo_grid;
  init_topology(grid, top, mpp, conf);

  /*param, using old charmm param for data structure compatliance*/
  if (rank == 0) puts("  params");
  charmm_param_t param;
  init_params(param, *grid.type_map, top, conf);

  /* velocity */
  create_velocity(&grid, conf.temperature, &mpp);
  int ntypes = grid.type_map->size();

  // if (rank == 0) {
  //   printf("Correcting replicated tags\n");
  // }
    
  // if (rank == 0)
  //   puts("sorting cells");

  if (rank == 0)
    puts("initialize comm buffer");
  comm_init_buf(&mpp, &grid);
  /*grid for MSM*/
  msm_grid_t mgrid;
  if (conf.coultype == COUL_MSM) {
    msm_grid_init(&mgrid, &mpp, &conf.msmlevels, conf.msmorder, grid.rcut, &grid.skin, conf.coulconst);
  }
  forward_comm_most(&grid, &mpp);

  MPI_Barrier(MPI_COMM_WORLD);

#ifdef __sw__
  if (rank == 0)
    puts("initialize arch dependent data");
  swinit(&grid);
#endif
  signal(SIGSEGV, sigsegv);
  /**/
  if (rank == 0) {
    printf("Initialize bond in cells\n");
  }

  //find_bonds(&grid);

  ff_t ff;
  ff.mpp = &mpp;
  ff.use_fdotr = M_NONB | M_BOND;
  ff.present = M_NONB | M_BOND; // | M_LONG;
  ff.grid = &grid;
  ff.longgrid = &mgrid;
  ff.param = &param;
  ff.bonded = compute_listed_forces;
  ff.nonbonded = (nonbonded_force_t*)charmm_nonbonded;
  // auto spring = new spring_com;
  // spring->mol = 0;
  // spring->K = 0.00625;
  // ff.enhs.append(spring);
  if (conf.coultype == COUL_MSM) {
    ff.present |= M_LONG;
    ff.longrange = (longrange_force_t*)msm_compute;
  }

  // rigid_type_table rigid_table(elemtab.cnt);

  grid.topo->update(&mpp, &grid);
  real dtfsq = conf.dt * conf.dt * 0.5 / 48.88821291 / 48.88821291;

  rigid_t rigid;
  if (conf.rigidtype == RIGID_NONE) {
    rigid.mask = 0;
    rigid.setup = 0;
  } else if (conf.rigidtype == RIGID_SHAKE) {
    rigid.mask = RIGID_POST_FORCE;
    rigid.setup = shake_setup;
    rigid.post_force = shake_post_force;
  }
  else if (conf.rigidtype == RIGID_LINCS) {
    rigid.mask = RIGID_POST_FORCE;
    rigid.setup = lincs_setup;
    rigid.post_force = lincs_post_force;
  }
  if (conf.tempcouple == TC_NH) {
    ff.temp_integrate = (temp_integrate_t*)nh_temp_integrate;
    nh_tstat_t tstat;
    nh_init(&tstat, &grid, &conf);
    ff.tstat = &tstat;
    ff.present |= M_TEMP_INTEGRATE;
  }
  param.nonb.order = conf.msmorder;
  integrate_verlet(0, &ff, &rigid, &mpp, &conf);
#ifdef __sw__
  swfinal();
#endif
  MPI_Finalize();
  char fn[PATH_MAX];
  sprintf(fn, "perf.%06d.log", rank);
  FILE *perff = fopen(fn, "w");

  timer_fprint(perff);
  fclose(perff);
  if (rank == 0) {
    stop_prof();
    timer_fprint(stdout);
    esmd::print_memusage();
  }
  return 0;
}
